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We study the quantum tunneling dynamics of many-body entangled solitons composed of ultracold 
bosonic gases in ID optical lattices. A bright soliton, confined by a potential barrier, is allowed 
to tunnel out of confinement by reducing the barrier width and for varying strengths of attractive 
interactions. Simulation of the Bose Hubbard Hamiltonian is performed with time-evolving block 
decimation. We find the characteristic 1/e time for the escape of the soliton, substantially different 
from the mean field prediction, and address how many-body effects like quantum fluctuations, 
entanglement, and nonlocal correlations affect macroscopic quantum tunneling; number fluctuations 
and second order correlations are suggested as experimental signatures. We find that while the 
escape time scales exponentially in the interactions, the time at which both the von Neumann 
entanglement entropy and the slope of number fluctuations is maximized scale only linearly. 



Tunneling is one of the most pervasive concepts in 
quantum mechanics and is essential to contexts as diverse 
as biophysics [l|, the a-decay of nuclei, vacuum states in 
quantum cosmology [2[ and chromodynamics [3[ . Macro- 
scopic quantum tunneling (MQT), the aggregate tunnel- 
ing behavior of a quantum many-body wavefunction, has 
been demonstrated in many condensed matter systems [4| 
and is one of the remarkable features of Bose-Einstein 
Condensates (BECs). For example, predictions for MQT 
in BECs range from cold atom Josephson rings [5[ to col- 
lapsing BECs [6[_, and MQT has been observed in double 
well potentials [7]. MQT has mainly been treated under 
semiclassical approximations such as JWKB and instan- 
ton methods, while more recently significant progress has 
been made towards a more general many-body picture 
via multi-configurational Hartree-Fock theory [8]. We 
present the first fully many-body entangled dynamical 
study of the quantum tunneling escape problem. 

In this Letter, we study how an initially trapped bright 
soliton composed of an ultracold bosonic gas in an optical 
lattice escapes via quantum tunneling through a barrier. 
In particular, we find that the maximum of both the 
von Neumann entanglement entropy and the slope of the 
number fluctuations between the initially trapped region 
behind the barrier and the escape region increases lin- 
early with the interaction strength, while the escape time 
tesc, i-e., the time at which the average number of remain- 
ing particles falls to 1/e of its initial value, increases expo- 
nentially with interactions. The optical lattice is required 
to use time-evolving block decimation (TEBD) [9], and 
also serves to increase quantum depletion and drive the 
system towards the strongly non-semi-classical regime. 

In the mean field limit, the macroscopic wavefunc- 
tion of the condensate is well-described by the Gross- 
Pitaevskii or nonlinear Schrodinger equation (NLS). The 
quantum many-body escape problem has already been 
studied in this context, where it has been found to have 
quite different features from single-particle quantum tun- 
neling, including a tunneling time which is not simply 
the inverse of the JWKB tunneling rate [loj and a non- 



smooth dynamical behavior referred to as "blips," in 
which particles escape through the barrier in bursts [ll[ • 

The NLS admits bright and dark solitonic solutions; 
we focus on the case of attractive interactions, of strong 
interest in recent experiments due to the newly demon- 
strated ability to change the interaction sign and strength 
over seven orders of magnitude [12). Solitons are local- 
ized, robust waves that propagate over long distances 
without changing shape, due to the nonlinear ity in the 
NLS, which counteracts dispersion. However, away from 
the mean field limit, quantum fluctuations affect the dy- 
namics and stability of solitons [1311. Bright solitons in 
BECs have already been observed [l4j , including in opti- 
cal lattices [l5| , so that our predictions and ideas can be 
tested with present experimental apparatus. In particu- 
lar, the small number of atoms we work with, from a few 
to 70, can be created in a 2D array of ID systems [l6j |. 
Proposed applications of bright solitons include a pulsed 
atomic laser [l7j. In such applications, MQT, quantum 
entanglement, and quantum fluctuations are all critical 
contributors to the overall dynamics and stability of soli- 
tons, and must be taken into account. 

Consider a system of N bosons at zero temperature in 
the canonical ensemble held in a ID homogeneous lat- 
tice of L sites, with box boundary conditions, similar to 
Ref. [16]. The lattice is sufficiently deep to allow us to 
evoke a lowest Bloch band tight-binding approximation. 
Then the system is well-described by a discretized ver- 
sion of the full many-body Hamiltonian, i.e., the Bose 
Hubbard Hamiltonian (BHH), 
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In Eq. (fT|), J is the energy of hopping and U < deter- 
mines the on-site two-particle interactions. An external 
rectangular potential barrier, of width w and height h, is 
given by V" ext . The field operator b\ (hi) creates (annihi- 
lates) a boson at the zth site and hi = b\bi. We will work 
in hopping units: J = 1 and time t in units of H/J. 
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FIG. 1: Initial State. Bright soliton formed by relaxation in 
imaginary time with a barrier of height h and initial width wi 
(dashed line). Before real time propagation the barrier is re- 
duced to width w (solid red line) so the soliton can commence 
macroscopic quantum tunneling. 
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To evolve the BHH in real and imaginary time we use 
TEBD TEBD is a matrix product state numeri- 
cal method that time evolves Eq. (pQ) on an adaptively 
reduced Hilbert space, given that the system is lowly 
entangled. TEBD is necessary because MQT is trig- 
gered by long range density fluctuations, and thus poorly 
modeled by mean field theory [l8[. Instanton meth- 
ods offer another approach towards calculating tunneling 
rates within a semiclassical approximation [19], but are 
rendered inaccurate for larger interaction strengths [20| , 
whereas TEBD suffers from no such limitations. 

To describe the system from a mean field perspective, 
the discrete NLS (DNLS) may either be obtained via dis- 
cretization of the NLS or from a mean field approxima- 
tion of the BHH. In the latter case, one can propagate 
the field operator hi forward in time using the BHH in 
the Heisenberg picture: ihd t bi = [bi,H]. Assuming the 
many-body state is a product of Glauber coherent states, 
(b\bibi) ^*^t^t, where ^ (U), leads to the DNLS: 

ihiia = -Jtyi+i + tl>i-i) + g\1>i\ 2 il>i + Vr^i- (2) 

In Eq. (j2j), the condensate order parameter, is nor- 
malized to the number of particles, N — Y^f=i IV^I 2 - 
Mean field simulations are performed using a pseudospec- 
tral, fourth-order Runge-Kutta, adaptation of Eq. (j2j). 
The BHH approaches the DNLS in the mean field limit 
N — > oo, U —^0, NU/J = const. We emphasize that 
both the BHH and the DNLS are single band models, 
valid when the soliton covers many sites; a continuum 
limit is possible for vJ = const., v = N/L — »• and 
J — >• oo; however, a continuum limit would restrict us 
numerically to very small numbers of particles (2l| . 

We initialize the many-body wavefunction via imagi- 
nary time relaxation as a bright soliton trapped behind 
the barrier as illustrated in Fig. HJ We set F ext to height 
h = 0.05 and width wi, effectively reducing the system 
size. At t = 0, in real time, the barrier is decreased to 
width w, where w is typically one to five sites, such that 
the soliton can escape on a time scale within reach of 
TEBD simulations. Attractive interactions U < < h 



FIG. 2: Many-Body Tunneling, and Calculation of Decay 
Time. Average particle number per site (top row) and number 
of trapped (blue) and escaped (red) particles (bottom row), 
for barrier widths w = 1 (left), w = 2 (middle), and w = 4 
(right); the 1/e decay time is t™f = 96.9, 123.9, and 219.4, 
all ±1.25, respectively. 



ensure that tunneling is always quantum, not classical. 
We choose L large enough so that reflections from the 
box boundary at the far right do not return to the bar- 
rier in t < t esc . Evolving in real time, we first make a 
coarse observation of the dynamics of MQT in Fig. [2] by 
plotting the average particle number in different regions, 
in order to determine t esc . 

How do many-body predictions compare to mean field 
ones? We define and t^f as the mean field and 
many-body escape times, respectively. For fixed NU/J, 
w, and h, the DNLS gives the same result independent of 
TV and U ; t^B _^ t MF only for N ^ and jj q. and 

w 2 h determines the barrier area. Figure [3] illustrates our 
exploration of this parameter space. The dynamics of 
MQT predicted by the DNLS and BHH differ strongly. 
Generally, the DNLS grossly under predicts t esc when 
N\U\/ J is sufficiently large. For example, in Fig. [3fc) 
for N\U\/J = 0.15, except for N = 1, the BHH predicts 
an increase in t esc over the DNLS, approaching a nearly 
constant value for v ~ 1/2 to 4. Although out of range of 
our simulations, we expect to subsequently decrease 
back down to tj^f , when v ^> 1. Escape times follow the 
same qualitative pattern for higher N\U\/J. The increase 
is due to number fluctuations, which are not permitted 
by mean field theory; particles spend more time further 
from the barrier at the soliton peak, and thus take longer 
to tunnel. 

Systematic error in results from the Schmidt trun- 
cation used in TEBD [9], %, the truncation in the on- 
site Hilbert space dimension, and the time resolution 
at which we write out data, St. The hardest many- 
body measures to converge, such as the block entropy, 
at x — 35 have an error < 10 -3 for N = 70, and were 
checked up through x — 55; due to small U and effec- 
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FIG. 3: Many Body vs. Mean Field Escape Time Predictions. 
(a)-(b) Dependence of £^<? on barrier area and particle num- 
ber for (a) NU/J = -0.15 and (b) -0.3. (c) t™£ plateaus 
for 10 to 70 particles as shown for NU/J = —0.15. (d) In the 
plateau region of N — 40, £^f<? significantly differs from 
for a range of barrier areas and interaction strengths. Curves 
are a guide to the eye, points represent actual data. 
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FIG. 4: Many-body Quantum Measures, (a) Average number 
at the peak of the soliton shows bursts of particles [ill ], (b) 
Fluctuations in the number of trapped particles increases with 
\U\. (c) Universal curve for the entropy of entanglement vs. 
the average number of trapped particles, (d) Exponential vs. 
linear increase in escape vs. many body times as a function of 
interactions. (Inset) Block entropy and dfi/dt closely follow 
each other, here for \U\ = 0.06. The key applies to panels 
(a)-(c) and all plots treat N = 6. 



tive system size, much lower x 1S required than usual 
in TEBD. Error bars in Fig. [3fd) are due solely to St] 
there is additional error from our Suzuki- Trotter expan- 
sion which is smaller than ^-induced error. For up to 
N = 20 we have not truncated <i, but for larger N up to 
70, we truncated to d = 18. A lower truncation results 
in decreased t^f , e.g. by 10% for d = 5, NU/J = -0.1, 
and TV = 10, even though max((n)) < 1, since more 
weight is given to spread-out Fock states. The attractive 
BHH requires much higher d than the repulsive BHH, 
since U < increases number fluctuations in high den- 
sity regions, i.e., at the soliton peak. The BHH also has 
a number of sources of systematic error, the most impor- 
tant of which is virtual fluctuations to the second band; 
however, since we compare single-band DNLS to single- 
band BHH this does not effect our comparison. 

In Fig. (D[a) we plot the average number at the peak 
of the soliton. There are points in time when the num- 
ber density exhibits steep exponential decay, and others 
during which it is nearly constant, similar to the den- 
sity bursts found by Dekel et al. [Ill; thus their predic- 
tions are correct even in the many-body regime. The 
first burst is independent of U. The initial flat horizon- 
tal region, at t c± 25, originates from initially left moving 
particles that are reflected off the leftmost infinite bound- 
ary and return back to the barrier. All subsequent de- 
viations from exponential decay appear to be dependent 
on U. In the coarser measure we find exponential 



scaling dominates for stronger interactions, as shown in 
Fig. [3(d). The dependence is exponential for two reasons: 
the many-body wavefunction tends to have large number 
fluctuations at the soliton peak, keeping particles away 
from the barrier; and the averaged density creates an ef- 
fective potential which increases the effective barrier size, 
as in mean field theory. Results in Fig. [4] are for N = 6; 
we found qualitatively similar results for up to N = 20, 
although simulations are limited in the large \U\ regime. 

To characterize the quantum nature of MQT, in 
Fig. [4(b) we plot the fluctuations in the number of parti- 
cles behind the barrier, f t = ((Af ) - (7V Z ) 2 )/ (TV/), where 
Ni is the number of particles to the left of site I and / 
is taken at the outer edge of the barrier. Once MQT 
commences, the maximum value of f\ in time increases 
with \U\ because number densities just outside the bar- 
rier have more influence to "pull" additional particles 
through the barrier, and vice versa. For early times, less 
than t = 50 in Fig. [4(b), // increases faster for smaller 
values of \U\ because the initial soliton is wider, but in- 
teractions take over shortly thereafter. 

Song et al. [22[ have shown that in ID conformal 
systems for which there is a conserved quantity, such 
as particle number, the variance of the fluctuations in 
that quantity between two subsystems A, B scales with 
the von Neumann entanglement entropy between A, B. 
Of particular interest to MQT is the von Neumann 
block entropy characterizing entanglement between the 
remaining particles and the escaped particles, Si = 
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FIG. 5: Time- dependence of Density- density Correlations. 
(a)-(c) shows correlations between trapped and escaped 
particles. The barrier, indicated by dashed lines, breaks up 
negatively-correlated regions (red); shown are time slices at 
t = 0, 62, 125. (d) Quantum depletion also grows rapidly. 



— Ti(pi log pi), where pi is the reduced density matrix for 
the well plus barrier. The key features of Si are illus- 
trated in a universal curve in Fig.|4^c): on the lower right 
side tunneling has not yet commenced; Si maximizes part 
way through the tunneling process in the center of the 
curve; and Si then decreases again to the left as the parti- 
cles finish tunneling out. Defining t s as the time at which 
Si is maximized and tf as the time at which the slope of 
the number fluctuations, dfi/dt, is maximized, we find t s 
and tf both increase linearly with U, in contrast to t^p, 
as shown in Fig. HJd). Moreover, Si and dfi/dt follow 
the same general trends in time, as shown in the inset 
of Fig. H(d). The two do not scale precisely, as dfi/dt is 
distorted by the density bursts illustrated in Fig. HJ^a). 

Another experimental signature is the density-density 

(2) 

correlations, g\- } = (fiihj) — (fii)(nj), extractable from 
noise measurements [23|; is zero in mean field the- 
ory. As customary, we subtract off the large diago- 
nal matrix elements of g^ to view the underlying off- 
diagonal structure. In Fig. Efa)-(c) we show for 
N = 40, NU/J = -0.015, and w = 2, dividing up the 
system to observe correlations between the three phys- 
ical regions: trapped, under the barrier, and escaped. 
We initially observe near-zero correlations everywhere 
except near the soliton peak. At t s , g^ shows many 
negatively-correlated regions (g^ < 0) which are bro- 
ken up by the potential barrier. In Fig. [5jd) we also 
show, for N — 20, rapidly growing quantum depletion 
D = 1 - (Em=2 A m)/(Em=i X m) where {A m } are the 
eigenvalues of the single particle density matrix (b\bj). 
This growth in D emphasizes the many-body nature of 
the escape process. 

In conclusion, we have performed quantum many-body 



simulations of the macroscopic quantum tunneling of 
bright solitons using TEBD to time-evolve the attrac- 
tive Bose-Hubbard Hamiltonian. We found strong devi- 
ations from mean field predictions. The escape time was 
shown to increase exponentially with interactions while 
both block entropy and the slope of number fluctuations 
maximized at a time which scaled linearly; entropy gen- 
erally followed closely the slope of number fluctuations, 
suggesting a dynamical extension of the static concepts 
of Song et al 22]. Our study shows that many-body ef- 
fects in macroscopic quantum tunneling can be observed 
via number fluctuations and density-density correlations 
as well as the increased escape time. 
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